RESULTS excluding rare and separating species by groups

Running models with different datasets:

Using only the 5x5 m quadrat scale

#Figures colors/labels
cores <- c( quadrat =  "#c00000", 
           `sp:time`= "#203864",
            time = "#ed7d31",
           `quadrat:time` = "#843c0c", 
           `quadrat:sp`= "#b666d2", 
            sp = "#70ad47", 
            Residual="#ffc000")

labvit <- as_labeller(c(grow = "Growth", mort="Mortality", 
                        rec="Recruitment"))
grpTlab = as_labeller(c(quadrat = "space", `sp:time`="species x time", 
                        `quadrat:time`="space x time", `quadrat:sp`= "species x space",
                        sp="species", Residual = "residual",
                        time= "time"))
grpvitT <- as_labeller(c(grow = "Growth", mort="Mortality", 
                         rec="Recruitment",
                        quadrat = "space", `sp:time`="species*time", 
                         `quadrat:time`="space*time", `quadrat:sp`= "species*space",
                         sp="species", Residual = "residual",
                        time= "time"))
labplot = as_labeller(c(quadrat = "space", 
                        `sp:time`="species x time",
                        `quadrat:time`="space x time", 
                        `quadrat:sp`= "species x space",
                        sp="species", 
                        Residual = "residual",
                        time= "time",
                fus="Fushan", luq = "Luquillo", lam="Lambir", 
                pas = "Pasoh", bci="Barro Colorado"))

Data

Growth

grp <- c("all","exclude","regroup")
path1 = here("models_outputs/models_time/")
tudo <- list()
for (i in 1:length(grp)){
  path <- paste0(path1, grp[i], "/grow/table")
  files = list.files(path)
  # excluding wyw
  files = files[grep("wyw", files, invert = T)]
  all <-  lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
  names(all) <- substr(files,1,nchar(files)-6)
  tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")

#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$data,tes$fplot)/8
##          
##           bci fus lam luq pas
##   all      10  10  10  10  10
##   exclude  10  10  10  10  10
##   regroup  10  10  10  10  10
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"

grow <- tes %>% filter(term != "intercept") %>%
  left_join(intercepts[,c(1:4,9)], by= c("id", "sub","data", "fplot")) %>%
  group_by(id, sub, data, fplot) %>% 
  mutate(VPC = variance/sum(variance))  %>% ungroup() %>%
  mutate(term = fct_relevel(term,"time", "quadrat", "sp:time", "quadrat:time",
                            "quadrat:sp", "sp", "Residual")) %>%
  as.data.frame()

# divide the sd of the terms by the mean growth (intercept) for each dataset.
grow$stdev = grow$Estimate/grow$intercept

Mort

tudo <- list()
for (i in 1:length(grp)){
  path <- paste0(path1,grp[i], "/mort/table")
  files = list.files(path)  # excluding wyw
  files = files[grep("wyw", files, invert = T)]
  all <-  lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
  names(all) <- substr(files,1,nchar(files)-6)
  tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")

#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$data,tes$fplot)/8
##          
##           bci fus lam luq pas
##   all      10  10  10  10  10
##   exclude  10  10  10  10  10
##   regroup  10  10  10  10  10
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"

mort <- tes %>% filter(term != "intercept") %>%
  left_join(intercepts[,c(1:4,9)], by= c("id", "sub","data", "fplot")) %>%
  group_by(id, sub, data, fplot) %>% 
  mutate(VPC = variance/sum(variance))  %>% ungroup() %>%
  mutate(term = fct_relevel(term,"time", "quadrat", "sp:time", "quadrat:time",
                            "quadrat:sp", "sp", "Residual")) %>%
  as.data.frame()
mort$stdev <- mort$Estimate

Rec

tudo <- list()
for (i in 1:length(grp)){
  path <- paste0(path1, grp[i], "/rec/table")
  files = list.files(path)  # excluding wyw
  files = files[grep("wyw", files, invert = T)]
  all <-  lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
  names(all) <- substr(files,1,nchar(files)-6)
  tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")

#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$data,tes$fplot)/8
##          
##           bci fus lam luq pas
##   all      10  10  10  10  10
##   exclude  10  10  10  10  10
##   regroup  10  10  10  10  10
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"

rec <- tes %>% filter(term != "intercept") %>%
  left_join(intercepts[,c(1:4,9)], by= c("id", "sub","data", "fplot")) %>%
  group_by(id, sub, data, fplot) %>% 
  mutate(VPC = variance/sum(variance))  %>% ungroup() %>%
  mutate(term = fct_relevel(term,"time", "quadrat", "sp:time", "quadrat:time",
                            "quadrat:sp", "sp", "Residual")) %>%
  as.data.frame()
rec$stdev <- rec$Estimate

Combining all results and taking the mean stdev and VPC among subsamples

NOT taking into account the same dataset (sub)

comb <- bind_rows(list(grow=grow, mort=mort, rec=rec), .id="vital") %>%
  mutate(term = fct_relevel(term, "quadrat", "sp:time","time", "quadrat:time", 
                            "quadrat:sp", "sp","Residual")) %>%
  ungroup() %>%
  mutate(data = fct_relevel(data, "all","regroup", "exclude"))
mcomb <- comb %>% group_by(vital,data,fplot,q_size,term) %>%
  summarise(stdev = mean(stdev),
            VPC = mean(VPC),
            richness = round(mean(richness)))  %>%
  filter(q_size == "quad_5")

number of fplots with results

table(mcomb$vital, mcomb$data)/7
##       
##        all regroup exclude
##   grow   5       5       5
##   mort   5       5       5
##   rec    5       5       5

INCLUDING THE % of rare species

load(here("data/rare_species.Rdata"))

comb <- comb %>% left_join(rare, by="fplot")
mcomb <- mcomb %>% left_join(rare, by="fplot")

Comparing groups

mcomb %>%
  ggplot(aes(x=data, y=VPC)) +geom_col(aes(fill=term, color=term), alpha=0.7) +
  scale_fill_manual(values=cores) +
  scale_color_manual(values=cores) +
  facet_grid(vital~fplot) +
  theme(axis.text.x = element_text(angle=45, hjust=1))

fig <- comb %>% filter(vital=="grow") %>%
  ggplot(aes(x=data, y=VPC, col=term)) + 
  geom_quasirandom(alpha=0.5, groupOnX = TRUE) +
  facet_grid(fplot~term, labeller=labplot) +
  stat_summary(fun=median, size=0.3, fatten=1,col="black",
                         geom="crossbar", width=0.7)+
  scale_color_manual(values = cores) +
    xlab("") +
  theme(legend.position = "none",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45,hjust=1)) +
  labs(title= "VPC per subsample of 5ha - GROWTH")
fig

fig <- comb %>% filter(vital=="mort") %>%
  ggplot(aes(x=data, y=VPC, col=term)) + 
  geom_quasirandom(alpha=0.5, groupOnX = TRUE) +
  facet_grid(fplot~term, labeller=labplot) +
  stat_summary(fun=median, size=0.3, fatten=1,col="black",
                         geom="crossbar", width=0.7)+
  scale_color_manual(values = cores) +
    xlab("") +
  theme(legend.position = "none",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45,hjust=1)) +
  labs(title= "VPC per subsample of 5ha - MORT")
fig

fig <- comb %>% filter(vital=="rec") %>%
  ggplot(aes(x=data, y=VPC, col=term)) + 
  geom_quasirandom(alpha=0.5, groupOnX = TRUE) +
  facet_grid(fplot~term, labeller=labplot) +
  stat_summary(fun=median, size=0.3, fatten=1,col="black",
                         geom="crossbar", width=0.7)+
  scale_color_manual(values = cores) +
    xlab("") +
  theme(legend.position = "none",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45,hjust=1)) +
  labs(title= "VPC per subsample of 5ha - REC")
fig

mean over fplots

mcomb %>% group_by(vital,data,term) %>% summarise(VPC=mean(VPC)) %>%
  ggplot(aes(x=data, y=VPC)) +geom_col(aes(fill=term, color=term), alpha=0.7) +
  scale_fill_manual(values=cores) +
  scale_color_manual(values=cores) +
  facet_grid(~vital)+
  theme(axis.text.x = element_text(angle=45, hjust=1))

comb %>%
  ggplot(aes(x=data,y=VPC, col=data)) + geom_boxplot()+
  facet_grid(vital~term) +
  theme(axis.text.x = element_text(angle=45,hjust=1),
        legend.position = "none")

mcomb %>% mutate(xis = as.numeric(data)) %>%
  ggplot(aes(x=xis,y=stdev, col=fplot)) + geom_line() +
  facet_grid(vital~term, scales = "free") + 
  scale_x_continuous(breaks = 1:3, labels=levels(mcomb$data),
                     name="") +
  theme(axis.text.x = element_text(angle=45, hjust=1))

USING the same dataset for the comparisons

It means, for each forest, linking the same subsampled data used

fig <- comb %>% filter(vital=="grow") %>%
  ggplot(aes(x=as.numeric(as.factor(data)), y=VPC, col=sub)) + 
  geom_line() +
  facet_grid(fplot~term, labeller=labplot) +
  scale_x_continuous(breaks=1:3, labels=levels(comb$data))+
    xlab("") +
  theme(legend.position = "none",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45,hjust=1)) +
  labs(title= "VPC per subsample of 5ha - GROW")
fig

fig <- comb %>% filter(vital=="mort") %>%
  ggplot(aes(x=as.numeric(as.factor(data)), y=VPC, col=sub)) + 
  geom_line() +
  facet_grid(fplot~term, labeller=labplot) +
  scale_x_continuous(breaks=1:3, labels=levels(comb$data))+
    xlab("") +
  theme(legend.position = "none",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45,hjust=1)) +
  labs(title= "VPC per subsample of 5ha - mort")
fig

fig <- comb %>% filter(vital=="rec") %>%
  ggplot(aes(x=as.numeric(as.factor(data)), y=VPC, col=sub)) + 
  geom_line() +
  facet_grid(fplot~term, labeller=labplot) +
  scale_x_continuous(breaks=1:3, labels=levels(comb$data))+
    xlab("") +
  theme(legend.position = "none",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45,hjust=1)) +
  labs(title= "VPC per subsample of 5ha - rec")
fig

labplot2 = as_labeller(c(grow = "Growth", mort="Mortality", 
                        rec="Recruitment",
                fus="Fushan", luq = "Luquillo", lam="Lambir", 
                pas = "Pasoh", bci="Barro Colorado")) 
fg <- mcomb %>% filter(data %in% c("all", "exclude", "regroup")) %>% 
    mutate(data = fct_relevel(data, "all", "exclude", "regroup")) %>%
  ggplot(aes(x=data, y=VPC)) +
  geom_col(aes(fill=term, color=term), alpha=0.7) +
  scale_fill_manual(values=cores,labels=
                       c("Space","Species x time", "Time", "Space x time", "Species x Space","Species","Residual")) +
  scale_color_manual(values=cores, labels=
                       c("Space","Species x time", "Time", "Space x time", "Species x Space","Species","Residual")) +
  
  facet_grid(fplot~vital, labeller = labplot2) +
  theme(axis.text.x = element_text(angle=45, hjust=1))
fg

# png("figs/exclude_regroup_rare_VPC_time_bar2.png", height = 800, width=800, res=100)
# fg
# dev.off()
mcomb %>% filter(data %in% c("all", "exclude", "regroup")) %>% 
  mutate(x=as.factor(data)) %>% mutate(x= as.numeric(x)) %>%
ggplot(aes(x=x, y=stdev, col=fplot)) + geom_line() +
  facet_grid(vital~term, scales="free") +
  scale_x_continuous(breaks = 1:3, labels=c("all","regroup", "exclude"),
                     name="")+
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "right") +
  ggtitle("SD")

mcomb %>% filter(data %in% c("all","regroup", "exclude")) %>% 
  mutate(x=as.factor(data)) %>% mutate(x= as.numeric(x)) %>%
ggplot(aes(x=data, y=stdev, col=data)) + geom_boxplot() +
  facet_grid(vital~term, scales="free") +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none") +
  ggtitle("SD")

mcomb %>% filter(data %in% c("all","regroup", "exclude")) %>% 
  mutate(x=as.factor(data)) %>% mutate(x= as.numeric(x)) %>%
ggplot(aes(x=x, y=VPC, col=fplot)) + geom_line() +
  facet_grid(vital~term, scales="free") +
  scale_x_continuous(breaks = 1:3, labels=c("all","regroup", "exclude"),
                     name="")+
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "right") +
  ggtitle("VPC")

mcomb %>% filter(data %in% c("all","regroup", "exclude")) %>% 
  mutate(x=as.factor(data)) %>% mutate(x= as.numeric(x)) %>%
ggplot(aes(x=data, y=VPC, col=data)) + geom_boxplot() +
  facet_grid(vital~term, scales="free") +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none") +
  ggtitle("VPC")

Difference to all data

Absolute difference

dif <- mcomb %>% filter(data %in% c("all","regroup", "exclude")) %>%
  mutate(data = fct_relevel(data, "all", "exclude", "regroup")) %>%
  pivot_wider(id_cols = c(vital,fplot,term, p.rare),
                      names_from = data,
                      values_from = c(stdev,VPC, richness)) %>%
  mutate(stdev_difr.regroup = stdev_regroup - stdev_all,
         stdev_difr.exclude = stdev_exclude - stdev_all,
         VPC_difr.regroup   = VPC_regroup - VPC_all,
         VPC_difr.exclude   = VPC_exclude - VPC_all,
         prop.rare.sp = (richness_exclude/richness_all)) %>%
  select(1:4,14:18) %>%
  pivot_longer(5:8,  names_to = c(".value", "group"),
   names_pattern = "(.*)_(.*)") %>%
  mutate(group = substr(group,5,nchar(group)))
dif %>%
  ggplot(aes(x=group, y=stdev)) + 
  facet_grid(vital~term, scales= "free_x") +
  geom_boxplot() +
  geom_quasirandom(aes(col=fplot)) +
  ylab("SD group - SD all") + xlab("") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Absolute difference in SD") +
  theme(legend.position = "bottom",
         panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45, hjust=1))

dif %>%
  ggplot(aes(x=group, y=VPC)) + 
  facet_grid(vital~term, scales= "free_x") +
  geom_boxplot(outlier.alpha = 0) +
  geom_quasirandom(aes(col=fplot)) +
  ylab("VPC group - VPC all") + xlab("") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Absolute difference in VPC") +
  theme(legend.position = "bottom",
         panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45, hjust=1))

resdif <- dif %>% group_by(vital, term, group) %>% summarise(mean=mean(VPC, na.rm=T)) %>%
  pivot_wider(names_from = term, values_from = mean)

Relative differences in SD to all data

(group value - all data value)*100 / all data value

difr <- mcomb %>% filter(data %in% c("all","regroup", "exclude")) %>%
  mutate(data = fct_relevel(data, "all", "exclude", "regroup")) %>%
  pivot_wider(id_cols = c(vital,fplot,term, p.rare),
                      names_from = data,
                      values_from = c(stdev,VPC, richness)) %>%
  mutate(stdev_difr.regroup = (stdev_regroup - stdev_all)*100/stdev_all,
         stdev_difr.exclude = (stdev_exclude - stdev_all)*100/stdev_all,
         VPC_difr.regroup   = (VPC_regroup - VPC_all)*100/VPC_all,
         VPC_difr.exclude   = (VPC_exclude - VPC_all)*100/VPC_all,
         prop.rare.sp = (richness_exclude/richness_all)) %>%
  select(1:4,14:18) %>%
  pivot_longer(5:8,  names_to = c(".value", "group"),
   names_pattern = "(.*)_(.*)") %>%
  mutate(group = substr(group,5,nchar(group)))
difr %>%
  ggplot(aes(x=group, y=stdev)) + 
  facet_grid(vital~term, scales= "free_x") +
  geom_boxplot() +
  geom_quasirandom(aes(col=fplot)) +
  ylab("SD proportional diference (%)") + xlab("") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Relative difference in SD") +
  theme(legend.position = "bottom",
         panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45, hjust=1))

resdifr <- difr %>% group_by(vital, term, group) %>% summarise(mean=mean(stdev, na.rm=T)) %>%
  pivot_wider(names_from = term, values_from = mean)
resdifr %>%  kable(digits = 2)
vital group quadrat sp:time time quadrat:time quadrat:sp sp Residual
grow .exclude 11.65 -13.84 3.55 -3.95 1.86 -11.99 1.53
grow .regroup 13.69 -23.98 7.45 -6.76 0.31 -12.50 8.09
mort .exclude 5.47 -3.67 -0.48 0.30 -3.07 -5.19 0.00
mort .regroup 6.95 -5.41 0.07 -2.83 -0.94 -6.54 0.00
rec .exclude 21.03 -7.91 6.51 1.16 2.38 -7.74 0.00
rec .regroup 24.34 -10.11 5.37 1.79 -3.80 -10.61 0.00

Relative differences in SD

difr %>% group_by(vital, term, group) %>% summarise(mean=mean(stdev, na.rm=T)) %>%
  pivot_wider(names_from = c(vital,group), values_from = mean) %>%
  kable(digits=1)
term grow_.exclude grow_.regroup mort_.exclude mort_.regroup rec_.exclude rec_.regroup
quadrat 11.7 13.7 5.5 7.0 21.0 24.3
sp:time -13.8 -24.0 -3.7 -5.4 -7.9 -10.1
time 3.5 7.5 -0.5 0.1 6.5 5.4
quadrat:time -3.9 -6.8 0.3 -2.8 1.2 1.8
quadrat:sp 1.9 0.3 -3.1 -0.9 2.4 -3.8
sp -12.0 -12.5 -5.2 -6.5 -7.7 -10.6
Residual 1.5 8.1 0.0 0.0 0.0 0.0

Absolute differences in VPC

dif %>% group_by(vital, term, group) %>% summarise(mean=mean(VPC, na.rm=T)) %>%
  pivot_wider(names_from = c(vital,group), values_from = mean) %>%
  kable(digits=3)
term grow_.exclude grow_.regroup mort_.exclude mort_.regroup rec_.exclude rec_.regroup
quadrat 0.002 0.002 0.000 0.002 0.001 0.002
sp:time -0.007 -0.008 -0.001 -0.001 -0.005 -0.006
time 0.002 0.002 -0.001 0.000 0.009 0.011
quadrat:time 0.005 0.002 0.005 0.000 0.003 0.013
quadrat:sp 0.014 -0.006 -0.002 -0.001 0.007 -0.009
sp -0.042 -0.050 -0.020 -0.024 -0.021 -0.026
Residual 0.026 0.059 0.019 0.025 0.006 0.016

Comparing with Proportion of rare species

dif %>%
  ggplot(aes(x=p.rare, y=stdev, linetype=group)) + 
  geom_point(aes(col=fplot)) + 
  geom_smooth(method="lm", se=F) +
  facet_grid(vital~term) +
  ylab("Difference in SD to original data") +
  xlab("Prop of rare species") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Difference in SD") +
  theme(axis.text.x = element_text(angle=45, hjust=1))

difr %>%
  ggplot(aes(x=p.rare, y=stdev, linetype=group)) + 
  geom_point(aes(col=fplot)) + 
  geom_smooth(method="lm", se=F) +
  facet_grid(vital~term) +
  ylab("Proportional difference in SD to original data") +
  xlab("Proportion of rare species") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Relative difference in SD") +
  theme(axis.text.x = element_text(angle=45, hjust=1))

grpvitT <- as_labeller(c(grow = "Growth", mort="Mortality", 
                         rec="Recruitment",
                        quadrat = "Space", `sp:time`="Species x Time", 
                         `quadrat:time`="Space x Time", `quadrat:sp`= "Species x Space",
                         sp="Species", Residual = "Residual",
                        time= "Time"))
pm <- dif %>%
  ggplot(aes(x=p.rare, y=VPC)) + 
  geom_point(aes(col=fplot, shape=group), size=3) + 
  facet_grid(vital~term, labeller = grpvitT) +
  scale_x_continuous(breaks=c(50,55,60,65)) +
  scale_color_discrete(labels=c("Barro Colorado", "Fushan",
                                "Lambir", "Luquillo", "Pasoh"))+
  ylab("Difference in VPC to original data") +
  xlab("Proportion of rare species (%)") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  theme(legend.position="bottom",
         panel.spacing.y = unit(0.8, "cm"),
         panel.background = element_rect(fill="white", color="black"))+
  labs(tag="b)")
pm

# png("figs/FIG_S4.3b_exclude_regroup_rare_VPC_time.png", height = 600, width=1010, res=100)
# pm
# dev.off()
pmr <- difr %>%
  ggplot(aes(x=p.rare, y=stdev)) + 
  geom_point(aes(col=fplot, shape=group), size=3) + 
  facet_grid(vital~term, labeller = grpvitT) +
  scale_x_continuous(breaks=c(64,68,72)) +
  scale_color_discrete(labels=c("Barro Colorado", "Fushan",
                                "Lambir", "Luquillo", "Pasoh"))+
  ylab("Proportional difference in SD to original data") +
  xlab("Proportion of rare species (%)") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  theme(legend.position="bottom",
         panel.spacing.y = unit(0.8, "cm"),
          panel.background = element_rect(colour="lightgray"))
pmr

# png("figs/exclude_regroup_rare_SD_time_relativeDIFF.png", height = 600, width=1010, res=100)
# pmr
# dev.off()